勉強会の目的

学部生、修論生、博論生がRでよく使うデータ整形の手法について学ぶための勉強会です

Saplings@Ecologyは、学部生、修論生、博論生etc.向けの交流用Slackグループを設けています。登録はこちらから

テーマは「とりあえずRで何ができるかを知っておこう！」

知っている手法が複数あれば、選択の余地が生まれる
一度見ておけば、やり方がわからないときに「そういえば…」と思い出して突破口となりうる
「Rを絶対使え！」でも「使い方をマスターしろ！」でもなく、「こういう方法もあるよ」ということを伝えることが大事だと思います

統計処理ではなく、データ整形や集計、作図の方法についていくつか紹介します

目的（手段を増やす）に沿って研究で使いそうな処理をできるだけ使ってみます

※注意事項

この資料は「とりあえずTidyverseで出来ることを知る」という目的のもと、「使いそうだな」というTidyverseの関数をピックアップして紹介し、いくつかの実践的な例で処理を紹介する形式となっています。そのためTidyverseの設計思想や一連のパッケージの使い方を学ぶには不足があります。それらについてはR for Data Scienceのページで公式に詳しく解説されているので、そちらを見ていただければと思います。

またRそのものの使い方についてもかなり説明を省いている部分がありますので、それらについてはR-Tipsや関連書籍等で確認していただければと思います。

作図に関しても、ここでは「ggplot2を使ってとりあえず作図を行って、ある程度の体裁を整える」方法の解説を行っています。そのため、例としてして作成している図の体裁には、発表等に用いるにあたって不十分な部分があります。実際に作図を行われる際には雑誌の規定等を参照しながら各自で微調整を行ってください。

この資料で紹介する内容に基づいて起こる研究上の損害等については責任を負いかねますので、ここにあるコード等をご自身で用いられる場合には、公式のレファレンスでしっかり関数の仕様を理解する、自分自身でコードの確認を行うなどをした上で処理を実行してください。

勉強会の流れ

基礎編

Rとは？Rを使う理由・利点
データ整形の基礎
データ整形の実践
ggplot2を利用した作図の基礎

発展編

複数グループに対して同様の処理を行う（nest→mapによる反復処理）
ggplot2による地図の描画（sfパッケージを利用したGISデータの描画）

基礎編①Rとは？Rを使う理由・利点

Rとは？

R is a language and environment for statistical computing and graphics. It is a GNU project which is similar to the S language and environment which was developed at Bell Laboratories (formerly AT&T, now Lucent Technologies) by John Chambers and colleagues.

https://www.r-project.org/about.html　より

R core teamにより提供されているフリーの統計解析環境（ソフトウェア）
クロスプラットフォームで利用可能

特徴：

（まあまあわかりやすい）CUI操作
豊富なパッケージ
多様な形式の外部データの取り入れ、データの複数フォーマットでの出力に対応している

※Pythonとの違いは？　→個人的には、「Pythonはプログラミング言語、Rはコマンドで操作できるソフト」という印象を強く受ける。Rにおけるスクリプトは多くの場合、「プログラムを記述したもの」というよりも「データ処理の過程を記録したもの」である。

Rを使う上でとりあえず覚えておくといいこと

Rはコマンドで動きます。
基本的には、Rでの作業はオブジェクトとして格納されているデータを関数に渡して、何かしらの処理をし、出力を得ることの繰り返しです。
この「データ」には、CSVファイルなどから読み込んだ外部データやR上で生成したデータ（乱数など）などがあります。

R_flow

Rを使う理由・利点

R_merit

①繰り返し＆やり直しが楽

作業面

スクリプトをきちんと書いておけば、コードを実行するだけで処理を何回でも行うことができる
- 例えば、発表で用いるすべての図を一つのスクリプトを実行するだけで再出力することが可能
for文などを用いることにより容易に繰り返し実行が可能
- 例えば、「100通りのデータそれぞれに対して回帰分析を行う」という処理も簡単に実行できる
データが変わった場合にも同じコードで容易に再実行が可能
- 例えば追加調査などでデータが増えた場合にも、データの形式が変わっていなければ同じコードで処理を行うことができる

研究面

作業効率が上がる
試行錯誤をしやすくなることにより、研究手法をより最適化しやすくなる
- 例えば、複数の解析手法を試せるようになることで、より適した解析手法を選択できるようになる

②規格化を行いやすい

作業面

データの形式が同じでさえあれば、同じコードを違うデータに一律で適用し同じ形式の結果を出力することが可能
- 例えば、glm関数などはデータがデータフレーム形式で与えられれば、どんなデータに対しても処理を行うことができる

研究面

研究がどの程度一般化できるかを考える端緒となる
- 例えば関数定義などを考える過程で、「自分が扱っているデータは他の類似データとどこが違うのか」、「自分が行おうとしている処理は他のデータに対しても適用可能か」などについても考えることになる

③他人に伝えるのが簡単（≒再現しやすい）

作業面

どのようなデータ処理を行ったかは、コードを見れば一発でわかる
言葉や文で「こうこうこうやる」などと伝えるより、Rスクリプト一つを渡した方が圧倒的に楽
- 例えば、「どのような回帰分析を行ったか」の詳細を言葉のみで伝えることは非常に難しいが、処理そのものが記録されているスクリプトを共有すればそれを簡単に伝えることができる

研究面

発表の際にも方法の妥当性を保証する
- 例えば、他の人は公開されたコードを見てデータ処理の過程を確認し、実行することで処理を再現することができる
解析方法はRのコードを見ればすぐにわかる
- 例えば、「GLM解析を行った」の一言から尤度の算出方法や最尤推定値の算出方法等を知ることは出来ないが、model <- glm(data = data, ...)等のコードからはそれらの情報を含む、解析方法についての十分な情報を得ることができる

Tidyverseとは

The tidyverse is an opinionated collection of R packages designed for data science. All packages share an underlying design philosophy, grammar, and data structures.

https://www.tidyverse.org/　より

Tidyverseはデータを扱うために一貫性をもって作成されたパッケージ群
Tidyverseではデータフレームに対する操作を関数で行うことができる
- iris[, c(“Sepal.Length”, “Species”)] ⇔ select(iris, Sepal.Length, Species)
従来はデータフレームを入出力にできずに不便だった点も、Tidyverseでは解決されている

Ex. iris内における各種の観測数をカウント

table(iris$Species)

## 
##     setosa versicolor  virginica 
##         50         50         50

iris %>% count(Species)

##      Species  n
## 1     setosa 50
## 2 versicolor 50
## 3  virginica 50

Tidyverseを使うメリット

素のRを使うよりも圧倒的に楽＆わかりやすい（処理が明示的）

Ex.列を抜き出して、条件で選択

素のR

iris2 <- iris
iris2 <- iris2[, c("Sepal.Length", "Species")]
iris2 <- iris2[iris2[, "Sepal.Length"] < 5.0, ]

Tidyverse

iris2 <- iris %>% 
select(Sepal.Length, Species) %>% 
filter(Sepal.Length < 5.0)

生態学分野では多くの場合、調査で得たデータをデータフレームとしてRに読み込み、データ整形・解析・作図を行う。そのためデータフレーム志向で設計されたTidyverseは非常に使いやすく有用である（※個人的な意見です）

パイプについて

（簡単に言うと）処理を横につなげられるようにするツール

function(x, …) ⇔ x %>% function(…)

Ex. iris %>% select(Sepal.Length) ⇔ select(iris, Sepal.Length)

本当は段階的に別名オブジェクトに格納した方がコード的には確実だけど、 実際には連続で関数を適用したりするのでこっちの方が便利

unique(pull(filter(iris, Sepal.Length < 5.0), Sepal.Length))

## [1] 4.9 4.7 4.6 4.4 4.8 4.3 4.5

iris %>% filter(Sepal.Length < 5.0) %>% pull(Sepal.Length) %>% unique()

## [1] 4.9 4.7 4.6 4.4 4.8 4.3 4.5

パイプによってデータ処理のフローを、視覚的にもわかりやすく記述することができる

RStudioを使うメリット

RStudio

RStudio is an integrated development environment (IDE) for R. It includes a console, syntax-highlighting editor that supports direct code execution, as well as tools for plotting, history, debugging and workspace management.

https://www.rstudio.com/products/rstudio/#rstudio-desktop　より

RStudioとは

RStudioはR用の統合解析環境である。

RStudioとは

スクリプト↔︎コンソールの行き来が楽なので、スクリプトを書きやすい
コードに不慣れな人でも扱いやすい親切設計（色々なことがマウスクリックでもできる）なので、細かいことに煩わされずに済む
コードの色分けやオートフィル機能等もとても便利

基礎編②データ整形の基礎

準備

作業ディレクトリの設定

データや出力結果を置いておくためのディレクトリ（フォルダ）を指定する。

Rstudio上では Session > Set Working Directory > Choose Directory で設定できる。

今回はこちらのページでCode>Download ZIPと進み、ダウンロードして解凍したフォルダを作業ディレクトリにすることで以下コードを円滑に実行できる。

Tidyverseの読み込み

パッケージの読み込みはlibrary関数で行う。もしTidyverseがインストールされていない場合には、install.package("tidyverse")でパッケージをインストールする。

library(tidyverse)
# インストールは
# install.package("tidyverse")

使うデータ：島嶼における木本植物の種多様性に関するデータ

今回は以下のデータペーパーで提供されているデータを例として使用する。

Schrader, J., Moeljono, S., Tambing, J., Sattler, C., & Kreft, H. 2020. A new dataset on plant occurrences on small islands, including species abundances and functional traits across different spatial scales. Biodiversity Data Journal 8:.

URL: https://bdj.pensoft.net/article/55275/instance/5664007/

Use license: Creative Commons Public Domain Waiver (CC-Zero)

インドネシアのRaja Ampat諸島に位置する60の小規模な島における木本植物のデータセット。各島に2x2mのプロットからなるトランセクトを島の大きさに応じて設置し、プロット内に根を張る胸高直径2㎝以上のすべての種を記録したデータである。その他に、島についての詳細データや植物種の機能的特性のデータ、種名のリスト、GBIF上に記録されている種の出現記録のデータが付属している。

oo_412824

Schrader et al. (2020)　より

各データは以下のファイルに格納されている（元データを微調整し、CSV形式で保存したもの）：

Community_data.csv（プロット内の出現種のデータ）
Island_data.csv（島の詳細データ）
Plant_functional_trait_data.csv（植物の機能的特性データ）
Species_data.csv（種名リスト）
Plant_occeurences.csv（GBIF上の観察記録）

調査、データの詳細については詳しくはリンク先を参照のこと。

以下の項目名に示されている関数名にはレファレンスページへのリンクを貼ってあります。それぞれの関数を使う際には、ぜひ一度見てみてください（Rから呼び出せる関数ヘルプと同じ内容です）。

データの読み込み　read_csv

read_scv関数はread.csv関数とほぼ同等の関数であり、目的とするCSVファイル（やその他の形式のテキストファイル）のパスを渡すとデータの読み込みを行ってくれる。ただし、この際に読み込まれたデータはTibble形式で格納されることに注意。

※data.frameとtibbleは厳密には違いますが、この資料ではわかりやすさのために両方とも「データフレーム」と呼んで区別していません。

community <- read_csv("data/Community_data.csv")
island <- read_csv("data/Island_Data.csv")
trait <- read_csv("data/Plant_functional_trait_data.csv")
occurence <- read_csv("data/Plant_occurrences.csv")
species <- read_csv("data/Species_data.csv")

読み込んだデータの中身は以下のようになっている：

community

## # A tibble: 2,215 x 6
##    island_ID transect_ID plot_ID    species_ID DBH_cm tree_height_m
##    <chr>     <chr>       <chr>      <chr>       <dbl>         <dbl>
##  1 GB1       GB1_T2      GB1_T2_ST5 species_12      2           3  
##  2 GB1       GB1_T5      GB1_T5_ST2 species_25      2           3.1
##  3 GB3       GB3_T1      GB3_T1_ST3 species_22      2           2.8
##  4 GB3       GB3_T3      GB3_T3_ST3 species_3       2           2.5
##  5 GB3       GB3_T3      GB3_T3_ST5 species_8       2           2.2
##  6 GB3       GB3_T4      GB3_T4_ST5 species_5       2           3.1
##  7 GB5       GB5_T1      GB5_T1_ST1 species_35      2           2.2
##  8 GB6       GB6_T2      GB6_T2_ST5 species_38      2           2.4
##  9 GB7       GB7_T1      GB7_T1_ST1 species_8       2           3  
## 10 GB7       GB7_T1      GB7_T1_ST2 species_8       2           2.6
## # ... with 2,205 more rows

colnames(community)

## [1] "island_ID"     "transect_ID"   "plot_ID"       "species_ID"   
## [5] "DBH_cm"        "tree_height_m"

island

## # A tibble: 60 x 12
##    island_ID island_coordina~ island_area island_perimeter distance_Gam
##    <chr>     <chr>                  <dbl>            <dbl>        <dbl>
##  1 GB1       -0.520461, 130.~     4774.              303.          59.3
##  2 GB2       -0.517484, 130.~        7.29             10.2         56.6
##  3 GB3       -0.517587, 130.~     2330.              184.         172. 
##  4 GB4       -0.517875, 130.~        8.06             11.3        192. 
##  5 GB5       -0.517319, 130.~       20.3              16.6        136. 
##  6 GB6       -0.515251, 130.~      317.               85.1        382. 
##  7 GB7       -0.515023, 130.~     1575.              150.         345. 
##  8 GB8       -0.516626, 130.~     1264.              175.         107. 
##  9 GB9       -0.514338, 130.~     1716.              152.         400. 
## 10 GB10      -0.516784, 130.~      121.               46.4        115. 
## # ... with 50 more rows, and 7 more variables: buffer_area_1000m <dbl>,
## #   tree_basal_area <dbl>, species_number <dbl>, soil_depth_mean <dbl>,
## #   leaf_litter_cover <dbl>, no_transects <dbl>, no_plots <dbl>

colnames(island)

##  [1] "island_ID"          "island_coordinates" "island_area"       
##  [4] "island_perimeter"   "distance_Gam"       "buffer_area_1000m" 
##  [7] "tree_basal_area"    "species_number"     "soil_depth_mean"   
## [10] "leaf_litter_cover"  "no_transects"       "no_plots"

trait

## # A tibble: 57 x 13
##    species_ID chlorophyll_SPAD chlorophyll_mod fruit_mass seed_mass   LMA
##    <chr>                 <dbl>           <dbl>      <dbl>     <dbl> <dbl>
##  1 species_1              50              59.2       0.07     14.4   0.91
##  2 species_2              47.4            54.7       0.22     30.4   1.94
##  3 species_3              59.2            77.3       0.14     41.5   1.23
##  4 species_4              46.0            52.5       0.13      0.51  1.47
##  5 species_5              54.4            67.6      NA        14.7   1.17
##  6 species_6              43.4            48.2       0.51      2.51  1.43
##  7 species_7              41.4            45.1       0.01      6     0.56
##  8 species_8              57.9            74.6       0.59     88.6   1.72
##  9 species_9              56.8            72.3       0.2      74.5   0.79
## 10 species_10             NA              NA        NA        NA     1.17
## # ... with 47 more rows, and 7 more variables: leaf_area <dbl>,
## #   wood_density <dbl>, max_tree_height <dbl>, max_tree_height_3 <dbl>,
## #   leaf_N <dbl>, leaf_C <dbl>, leaf_P <dbl>

colnames(trait)

##  [1] "species_ID"        "chlorophyll_SPAD"  "chlorophyll_mod"  
##  [4] "fruit_mass"        "seed_mass"         "LMA"              
##  [7] "leaf_area"         "wood_density"      "max_tree_height"  
## [10] "max_tree_height_3" "leaf_N"            "leaf_C"           
## [13] "leaf_P"

occurence

## # A tibble: 373 x 17
##    id    basisOfRecord occurrenceID recordedBy eventDate islandGroup country
##    <chr> <chr>         <chr>        <chr>      <chr>     <chr>       <chr>  
##  1 plan~ LivingSpecim~ plant-occur~ Julian Sc~ 2016-06-~ Raja Ampat  Indone~
##  2 plan~ LivingSpecim~ plant-occur~ Julian Sc~ 2016-06-~ Raja Ampat  Indone~
##  3 plan~ LivingSpecim~ plant-occur~ Julian Sc~ 2016-06-~ Raja Ampat  Indone~
##  4 plan~ LivingSpecim~ plant-occur~ Julian Sc~ 2016-06-~ Raja Ampat  Indone~
##  5 plan~ LivingSpecim~ plant-occur~ Julian Sc~ 2016-06-~ Raja Ampat  Indone~
##  6 plan~ LivingSpecim~ plant-occur~ Julian Sc~ 2016-06-~ Raja Ampat  Indone~
##  7 plan~ LivingSpecim~ plant-occur~ Julian Sc~ 2016-06-~ Raja Ampat  Indone~
##  8 plan~ LivingSpecim~ plant-occur~ Julian Sc~ 2016-06-~ Raja Ampat  Indone~
##  9 plan~ LivingSpecim~ plant-occur~ Julian Sc~ 2016-06-~ Raja Ampat  Indone~
## 10 plan~ LivingSpecim~ plant-occur~ Julian Sc~ 2016-06-~ Raja Ampat  Indone~
## # ... with 363 more rows, and 10 more variables: countryCode <chr>,
## #   decimalLatitude <dbl>, decimalLongitude <dbl>, geodeticDatum <chr>,
## #   coordinateUncertaintyInMeters <dbl>, identificationQualifier <chr>,
## #   scientificName <chr>, kingdom <chr>, family <chr>, taxonRank <chr>

colnames(occurence)

##  [1] "id"                            "basisOfRecord"                
##  [3] "occurrenceID"                  "recordedBy"                   
##  [5] "eventDate"                     "islandGroup"                  
##  [7] "country"                       "countryCode"                  
##  [9] "decimalLatitude"               "decimalLongitude"             
## [11] "geodeticDatum"                 "coordinateUncertaintyInMeters"
## [13] "identificationQualifier"       "scientificName"               
## [15] "kingdom"                       "family"                       
## [17] "taxonRank"

species

## # A tibble: 57 x 5
##    species_ID Family     Species         Author           UNIPA_Voucher_ID      
##    <chr>      <chr>      <chr>           <chr>            <chr>                 
##  1 species_1  Rubiaceae  Ixora timorens~ Decne.           Schrader 12/13/82/99/~
##  2 species_2  Ebenaceae  Diospyros mari~ Blume            Schrader 14/30/40/80/~
##  3 species_3  Moraceae   Ficus microcar~ L.f.             Schrader 88/92/98/48/~
##  4 species_4  Moraceae   Ficus tinctoria G.Forst.         Schrader 11/97/155    
##  5 species_5  Primulace~ Myrsine rawace~ A. DC.           Schrader 35/95/173/101
##  6 species_6  Rubiaceae  Timonius sp. 1  <NA>             Schrader 20/151       
##  7 species_7  Euphorbia~ Macaranga dioi~ (G.Forst.) M・l.~ Schrader 4/176/112    
##  8 species_8  Myrtaceae  Eugenia reinwa~ (Blume) DC.      Schrader 10/121/159   
##  9 species_9  Micromelum Micromelum min~ (G.Forst.) Wigh~ Schrader 106/147      
## 10 species_10 Celastrac~ Pleurostylia o~ (Wall.) Alston   Schrader 58           
## # ... with 47 more rows

colnames(species)

## [1] "species_ID"       "Family"           "Species"          "Author"          
## [5] "UNIPA_Voucher_ID"

データ整形をする前に

Tidy_data

整然データについて詳しくはこちらのページを参照のこと。

Tidyverseによるデータ整形では、ほとんどの関数が列名を引数（関数に与える値）としてとる
つまり、「どの変数に対して何をしたいか」をはっきりさせておくことが大事
- 変数Aと変数Bの大小を比べたい！
- 変数Aを抜き出したい！
- 変数Aに関数を適用したい！

列の選択　select

selectはデータフレームから特定の列を抜き出すために用いる関数である。この際、抜き出された列がたとえ一列だったとしてもデータフレームとして出力される。ベクトルとして抜き出したいのであれば、後述するpullを用いる

community %>% 
select(island_ID, transect_ID, plot_ID, species_ID)

## # A tibble: 2,215 x 4
##    island_ID transect_ID plot_ID    species_ID
##    <chr>     <chr>       <chr>      <chr>     
##  1 GB1       GB1_T2      GB1_T2_ST5 species_12
##  2 GB1       GB1_T5      GB1_T5_ST2 species_25
##  3 GB3       GB3_T1      GB3_T1_ST3 species_22
##  4 GB3       GB3_T3      GB3_T3_ST3 species_3 
##  5 GB3       GB3_T3      GB3_T3_ST5 species_8 
##  6 GB3       GB3_T4      GB3_T4_ST5 species_5 
##  7 GB5       GB5_T1      GB5_T1_ST1 species_35
##  8 GB6       GB6_T2      GB6_T2_ST5 species_38
##  9 GB7       GB7_T1      GB7_T1_ST1 species_8 
## 10 GB7       GB7_T1      GB7_T1_ST2 species_8 
## # ... with 2,205 more rows

community %>% 
select(-DBH_cm, -tree_height_m)

## # A tibble: 2,215 x 4
##    island_ID transect_ID plot_ID    species_ID
##    <chr>     <chr>       <chr>      <chr>     
##  1 GB1       GB1_T2      GB1_T2_ST5 species_12
##  2 GB1       GB1_T5      GB1_T5_ST2 species_25
##  3 GB3       GB3_T1      GB3_T1_ST3 species_22
##  4 GB3       GB3_T3      GB3_T3_ST3 species_3 
##  5 GB3       GB3_T3      GB3_T3_ST5 species_8 
##  6 GB3       GB3_T4      GB3_T4_ST5 species_5 
##  7 GB5       GB5_T1      GB5_T1_ST1 species_35
##  8 GB6       GB6_T2      GB6_T2_ST5 species_38
##  9 GB7       GB7_T1      GB7_T1_ST1 species_8 
## 10 GB7       GB7_T1      GB7_T1_ST2 species_8 
## # ... with 2,205 more rows

条件による行の選択　filter

各列の値に基づく条件を用いて行を選択したい場合にはfilter関数を用いる。素のRで行う場合にはこの処理の記述は非常に煩雑になるため、filter関数は積極的に使っていきたい。

community %>% 
filter((island_ID == "GB7")&(transect_ID == "GB7_T1")&(plot_ID == "GB7_T1_ST2"))

## # A tibble: 4 x 6
##   island_ID transect_ID plot_ID    species_ID DBH_cm tree_height_m
##   <chr>     <chr>       <chr>      <chr>       <dbl>         <dbl>
## 1 GB7       GB7_T1      GB7_T1_ST2 species_8     2             2.6
## 2 GB7       GB7_T1      GB7_T1_ST2 species_8     2.3           2.5
## 3 GB7       GB7_T1      GB7_T1_ST2 species_1     2.3           2  
## 4 GB7       GB7_T1      GB7_T1_ST2 species_39    6.1           5.2

ベクトルの抜き出し　pull　新しい列に計算結果を格納　mutate　文字列に対する操作　str_関数

データフレーム中の列をベクトルとして抜き出したい場合にはpullを用いる。データフレームに、新たなデータを含む新たな列を追加したい場合にはmutateを用いる。多くの場合、既存の列の値をもとに新たな値を算出し、その値を格納する処理に用いられる。 str_関数はパッケージstringrに含まれる関数であり、ここで用いているstr_extractは与えられた文字列から特定パターン（正規表現も可能）と一致する部分を抜き出してくれる関数である。

occurence %>% pull(scientificName) -> occurence_scientificName
occurence_scientificName[1:10]

##  [1] "Glochidion castaneum Airy Shaw"      
##  [2] "Ixora timorensis Decne."             
##  [3] "Ficus pedunculosa Miq."              
##  [4] "Ficus prasinicarpa Elmer"            
##  [5] "Ficus microcarpa L.f."               
##  [6] "Ficus tinctoria G.Forst."            
##  [7] "Rapanea rawacensis (A. DC.) Mez"     
##  [8] "Aglaia elaeagnoidea (A.Juss.) Benth."
##  [9] "Severinia lauterbachii Swingle"      
## [10] "Planchonella obovata (R.Br.) Pierre"

occurence2 <- occurence %>% 
  mutate(species_name = str_extract(scientificName, pattern = "^[A-Za-z]+ [A-Za-z]+"))

occurence2 %>% pull(species_name) -> occurence2_scientificName
occurence2_scientificName[1:10]

##  [1] "Glochidion castaneum"   "Ixora timorensis"       "Ficus pedunculosa"     
##  [4] "Ficus prasinicarpa"     "Ficus microcarpa"       "Ficus tinctoria"       
##  [7] "Rapanea rawacensis"     "Aglaia elaeagnoidea"    "Severinia lauterbachii"
## [10] "Planchonella obovata"

列を分割　separate　日付の定義　lubridate::ymd

separateは特定列に含まれる文字列を、特定文字列を基準に分割し、2列に分けて格納する。ymdは日付を扱うためのパッケージlubridateに含まれる関数であり、「Year-Month-Day」形式で文字列として格納されている日付をDate型に変換してくれる。ここで用いている[パッケージ名]::[関数]という記法は特定のパッケージから1つだけ関数を呼び出したいときに便利。

occurence3 <- occurence2 %>% 
  separate(eventDate, sep = "/", into = c("Date1", "Date2")) %>% 
  mutate(
    Date1 = lubridate::ymd(Date1), 
    Date2 = lubridate::ymd(Date2)
  )

観測数のカウント　count

count関数は与えられた列内に含まれる各要素の観測数（＝データフレーム内での行数）を数え、データフレーム形式で出力してくれる。2つ以上の列名を与えた場合には、各要素の組み合わせそれぞれの観測数を数えてくれる。同様の処理はgroup_byとtallyを組み合わせることでも行うことができる。

community %>% 
  count(island_ID, species_ID)

## # A tibble: 390 x 3
##    island_ID species_ID     n
##    <chr>     <chr>      <int>
##  1 GB1       species_12     1
##  2 GB1       species_2      3
##  3 GB1       species_22     2
##  4 GB1       species_25     7
##  5 GB1       species_26     2
##  6 GB1       species_32     1
##  7 GB1       species_37     5
##  8 GB1       species_38     5
##  9 GB1       species_40     3
## 10 GB1       species_43     1
## # ... with 380 more rows

グループ化　group_by　と　集計　summerise

group_byは与えられた列の要素に基づいてデータフレーム内の行をグループ化する関数である。summariseは各グループに含まれるデータに対し関数を適用することで値を要約し、結果をデータフレームとして返す。

community %>% 
  group_by(species_ID) %>% 
  summarise(
    mean_DBH = mean(DBH_cm), 
    mean_height = mean(tree_height_m)
  )

## `summarise()` ungrouping output (override with `.groups` argument)

## # A tibble: 57 x 3
##    species_ID mean_DBH mean_height
##    <chr>         <dbl>       <dbl>
##  1 species_1      4.31        3.95
##  2 species_10     4.4         4.2 
##  3 species_11     2.86        2.88
##  4 species_12     2.86        3.28
##  5 species_13     2.69        3.46
##  6 species_15     4.9         4.52
##  7 species_16     3.7         3.6 
##  8 species_17     2.36        3.28
##  9 species_18     2.99        2.44
## 10 species_19     2.6         2   
## # ... with 47 more rows

データフレームの結合　full_join

full_joinはbyで与えられた列をキーとして2つのデータフレームを結合する関数である。full_joinは2つのデータフレームのキー列に含まれるすべての要素を含む（つまりデータフレームAとBをfull_joinで結合した場合、データフレームAに記録がない要素についての行ではデータフレームAに含まれていた列の値がすべてNAとなる）。

このように、full_joinによる結合ではデータフレームの行数が気付かぬうちに増えてしまう危険性がある。そのため、適宜左側結合のleft_join、右側結合のright_join等と使い分けるべきである。

species

## # A tibble: 57 x 5
##    species_ID Family     Species         Author           UNIPA_Voucher_ID      
##    <chr>      <chr>      <chr>           <chr>            <chr>                 
##  1 species_1  Rubiaceae  Ixora timorens~ Decne.           Schrader 12/13/82/99/~
##  2 species_2  Ebenaceae  Diospyros mari~ Blume            Schrader 14/30/40/80/~
##  3 species_3  Moraceae   Ficus microcar~ L.f.             Schrader 88/92/98/48/~
##  4 species_4  Moraceae   Ficus tinctoria G.Forst.         Schrader 11/97/155    
##  5 species_5  Primulace~ Myrsine rawace~ A. DC.           Schrader 35/95/173/101
##  6 species_6  Rubiaceae  Timonius sp. 1  <NA>             Schrader 20/151       
##  7 species_7  Euphorbia~ Macaranga dioi~ (G.Forst.) M・l.~ Schrader 4/176/112    
##  8 species_8  Myrtaceae  Eugenia reinwa~ (Blume) DC.      Schrader 10/121/159   
##  9 species_9  Micromelum Micromelum min~ (G.Forst.) Wigh~ Schrader 106/147      
## 10 species_10 Celastrac~ Pleurostylia o~ (Wall.) Alston   Schrader 58           
## # ... with 47 more rows

trait

## # A tibble: 57 x 13
##    species_ID chlorophyll_SPAD chlorophyll_mod fruit_mass seed_mass   LMA
##    <chr>                 <dbl>           <dbl>      <dbl>     <dbl> <dbl>
##  1 species_1              50              59.2       0.07     14.4   0.91
##  2 species_2              47.4            54.7       0.22     30.4   1.94
##  3 species_3              59.2            77.3       0.14     41.5   1.23
##  4 species_4              46.0            52.5       0.13      0.51  1.47
##  5 species_5              54.4            67.6      NA        14.7   1.17
##  6 species_6              43.4            48.2       0.51      2.51  1.43
##  7 species_7              41.4            45.1       0.01      6     0.56
##  8 species_8              57.9            74.6       0.59     88.6   1.72
##  9 species_9              56.8            72.3       0.2      74.5   0.79
## 10 species_10             NA              NA        NA        NA     1.17
## # ... with 47 more rows, and 7 more variables: leaf_area <dbl>,
## #   wood_density <dbl>, max_tree_height <dbl>, max_tree_height_3 <dbl>,
## #   leaf_N <dbl>, leaf_C <dbl>, leaf_P <dbl>

full_join(species, trait, by = c("species_ID" = "species_ID"))

## # A tibble: 57 x 17
##    species_ID Family Species Author UNIPA_Voucher_ID chlorophyll_SPAD
##    <chr>      <chr>  <chr>   <chr>  <chr>                       <dbl>
##  1 species_1  Rubia~ Ixora ~ Decne. Schrader 12/13/~             50  
##  2 species_2  Ebena~ Diospy~ Blume  Schrader 14/30/~             47.4
##  3 species_3  Morac~ Ficus ~ L.f.   Schrader 88/92/~             59.2
##  4 species_4  Morac~ Ficus ~ G.For~ Schrader 11/97/~             46.0
##  5 species_5  Primu~ Myrsin~ A. DC. Schrader 35/95/~             54.4
##  6 species_6  Rubia~ Timoni~ <NA>   Schrader 20/151              43.4
##  7 species_7  Eupho~ Macara~ (G.Fo~ Schrader 4/176/~             41.4
##  8 species_8  Myrta~ Eugeni~ (Blum~ Schrader 10/121~             57.9
##  9 species_9  Micro~ Microm~ (G.Fo~ Schrader 106/147             56.8
## 10 species_10 Celas~ Pleuro~ (Wall~ Schrader 58                  NA  
## # ... with 47 more rows, and 11 more variables: chlorophyll_mod <dbl>,
## #   fruit_mass <dbl>, seed_mass <dbl>, LMA <dbl>, leaf_area <dbl>,
## #   wood_density <dbl>, max_tree_height <dbl>, max_tree_height_3 <dbl>,
## #   leaf_N <dbl>, leaf_C <dbl>, leaf_P <dbl>

縦長変換　pivot_longer：列名と各列の値をそれぞれ一列に

pivot_longerはデータフレームを縦長変換するための関数であり、与えられた列の列名を一つの列に、各列の値をもう一つの列に格納する。特に読み込んだデータフレームが非Tidyな場合（例えば、各列に各月の値が入っている場合など）に、それらをTidyなデータに変換するために用いることができる。

island_pivot1 <- island %>% 
  select(island_ID, island_coordinates, island_area, island_perimeter)
island_pivot2 <- island_pivot1 %>% 
  pivot_longer(cols = c(island_area, island_perimeter), names_to = "Island_property", values_to = "Value")
island_pivot2

## # A tibble: 120 x 4
##    island_ID island_coordinates    Island_property    Value
##    <chr>     <chr>                 <chr>              <dbl>
##  1 GB1       -0.520461, 130.580800 island_area      4774.  
##  2 GB1       -0.520461, 130.580800 island_perimeter  303.  
##  3 GB2       -0.517484, 130.575020 island_area         7.29
##  4 GB2       -0.517484, 130.575020 island_perimeter   10.2 
##  5 GB3       -0.517587, 130.568534 island_area      2330.  
##  6 GB3       -0.517587, 130.568534 island_perimeter  184.  
##  7 GB4       -0.517875, 130.568515 island_area         8.06
##  8 GB4       -0.517875, 130.568515 island_perimeter   11.3 
##  9 GB5       -0.517319, 130.569511 island_area        20.3 
## 10 GB5       -0.517319, 130.569511 island_perimeter   16.6 
## # ... with 110 more rows

横長変換　pivot_wider：列内の要素を列名に

pivot_widerはデータフレームの横長変換を行うための関数であり、一つの列に含まれる要素を列名、もう一つの列に含まれる値を各列の値として引き出しデータフレームを横長に展開する。横長変換する際にはデータが存在しないセルも生じるが、その際にはそのセルの値はNAとなる。

pivot_longerの逆にあたる処理を行う関数。

island_pivot3 <- island_pivot2 %>% 
  pivot_wider(names_from = "Island_property", values_from = "Value")
island_pivot3

## # A tibble: 60 x 4
##    island_ID island_coordinates    island_area island_perimeter
##    <chr>     <chr>                       <dbl>            <dbl>
##  1 GB1       -0.520461, 130.580800     4774.              303. 
##  2 GB2       -0.517484, 130.575020        7.29             10.2
##  3 GB3       -0.517587, 130.568534     2330.              184. 
##  4 GB4       -0.517875, 130.568515        8.06             11.3
##  5 GB5       -0.517319, 130.569511       20.3              16.6
##  6 GB6       -0.515251, 130.568484      317.               85.1
##  7 GB7       -0.515023, 130.569425     1575.              150. 
##  8 GB8       -0.516626, 130.572132     1264.              175. 
##  9 GB9       -0.514338, 130.570002     1716.              152. 
## 10 GB10      -0.516784, 130.571887      121.               46.4
## # ... with 50 more rows

データフレームの書き出し　write_csv

write_csvはwrite.csvとほぼ同等の関数で、データフレームを指定されたファイル名のテキストファイル（CSV含む）に出力する。デフォルトのエンコーディングはUTF-8であるため、日本語が含まれる場合は注意。

write_csv(occurence3, file = "occurence3.csv")

基礎編③データ整形の実践

以下のデータを得る方法を考えてみよう！

ここまで見て来たTidyverseの関数を使って以下のような実践的な処理をどのように行うことができるか、考えてみよう。

紹介していない関数も使うので、適宜検索しながら進めよう。

種ごとのDBHの平均と標準誤差、種ごとの樹高の平均と標準誤差と種のデータ（species、trait）を結合したデータフレームを作成する
島ごとの総出現種数、総出現属数、総出現科数と島のデータを結合したデータフレームを作成する
調査プロット名×種名行列（＝PCAやDCAで使える形）を作成する

答えを見る前に考えてみてね

種ごとのDBHの平均と標準誤差、種ごとの樹高の平均と標準誤差と種のデータ（species、trait）を結合したデータフレームを作成する

group_byとsummariseを用いて、各種ごとに平均値に加えて標準誤差を算出する。得られたデータフレームをspecies_IDをキーとしてspecies、traitとfull_joinする。

community %>% 
  group_by(species_ID) %>% 
  summarise(
    mean_DBH = mean(DBH_cm), 
    se_DBH = sd(DBH_cm)/sqrt(length(DBH_cm)), 
    mean_height = mean(tree_height_m), 
    se_height = sd(tree_height_m)/sqrt(length(tree_height_m))
  ) -> species_dbh_height

## `summarise()` ungrouping output (override with `.groups` argument)

species_perperty <- full_join(species, trait, by = "species_ID") %>% full_join(species_dbh_height, by = "species_ID")

species_perperty

## # A tibble: 57 x 21
##    species_ID Family Species Author UNIPA_Voucher_ID chlorophyll_SPAD
##    <chr>      <chr>  <chr>   <chr>  <chr>                       <dbl>
##  1 species_1  Rubia~ Ixora ~ Decne. Schrader 12/13/~             50  
##  2 species_2  Ebena~ Diospy~ Blume  Schrader 14/30/~             47.4
##  3 species_3  Morac~ Ficus ~ L.f.   Schrader 88/92/~             59.2
##  4 species_4  Morac~ Ficus ~ G.For~ Schrader 11/97/~             46.0
##  5 species_5  Primu~ Myrsin~ A. DC. Schrader 35/95/~             54.4
##  6 species_6  Rubia~ Timoni~ <NA>   Schrader 20/151              43.4
##  7 species_7  Eupho~ Macara~ (G.Fo~ Schrader 4/176/~             41.4
##  8 species_8  Myrta~ Eugeni~ (Blum~ Schrader 10/121~             57.9
##  9 species_9  Micro~ Microm~ (G.Fo~ Schrader 106/147             56.8
## 10 species_10 Celas~ Pleuro~ (Wall~ Schrader 58                  NA  
## # ... with 47 more rows, and 15 more variables: chlorophyll_mod <dbl>,
## #   fruit_mass <dbl>, seed_mass <dbl>, LMA <dbl>, leaf_area <dbl>,
## #   wood_density <dbl>, max_tree_height <dbl>, max_tree_height_3 <dbl>,
## #   leaf_N <dbl>, leaf_C <dbl>, leaf_P <dbl>, mean_DBH <dbl>, se_DBH <dbl>,
## #   mean_height <dbl>, se_height <dbl>

島ごとの総出現種数、総出現属数、総出現科数と島のデータを結合したデータフレームを作成する

以下のようにselect、distinct、countを組み合わせることで、各島に存在する木本植物の種数、属数、科数を算出し、最後にisland_IDをキーに結合する。属数に関しては、str_extractで種名から属名を抽出し、それを参照してカウントする。

ただし、60島のデータを含むislandとそのうちの40島における調査データをもとに作成したデータフレームをfull_joinで結合しているため、island_propertyにはn_species等の列の値がNAになっている20行のデータが含まれることになる点には注意。

community_species <- left_join(community, species, by = "species_ID")

sp_per_island <- community_species %>% 
  select(island_ID, species_ID) %>% 
  distinct() %>% 
  count(island_ID, name = "n_species")
sp_per_island

## # A tibble: 40 x 2
##    island_ID n_species
##    <chr>         <int>
##  1 GB1              18
##  2 GB10              5
##  3 GB11              9
##  4 GB12             14
##  5 GB13             10
##  6 GB14             11
##  7 GB15              9
##  8 GB16              6
##  9 GB17              3
## 10 GB18              8
## # ... with 30 more rows

genera_per_island <- community_species %>% 
  mutate(Genera = str_extract(Species, pattern = "[A-Za-z]+")) %>% 
  select(island_ID, Genera) %>% 
  distinct() %>% 
  count(island_ID, name = "n_genera")
genera_per_island

## # A tibble: 40 x 2
##    island_ID n_genera
##    <chr>        <int>
##  1 GB1             17
##  2 GB10             5
##  3 GB11             9
##  4 GB12            13
##  5 GB13            10
##  6 GB14            11
##  7 GB15             9
##  8 GB16             6
##  9 GB17             3
## 10 GB18             7
## # ... with 30 more rows

family_per_island <- community_species %>% 
  count(island_ID, Family) %>% 
  select(island_ID, Family) %>% 
  distinct() %>% 
  count(island_ID, name = "n_family")
family_per_island

## # A tibble: 40 x 2
##    island_ID n_family
##    <chr>        <int>
##  1 GB1             14
##  2 GB10             5
##  3 GB11             6
##  4 GB12             9
##  5 GB13             8
##  6 GB14            10
##  7 GB15             8
##  8 GB16             6
##  9 GB17             2
## 10 GB18             7
## # ... with 30 more rows

island_property <- full_join(island, sp_per_island, by = "island_ID") %>% 
  full_join(genera_per_island, by = "island_ID") %>% 
  full_join(family_per_island, by = "island_ID")
island_property

## # A tibble: 60 x 15
##    island_ID island_coordina~ island_area island_perimeter distance_Gam
##    <chr>     <chr>                  <dbl>            <dbl>        <dbl>
##  1 GB1       -0.520461, 130.~     4774.              303.          59.3
##  2 GB2       -0.517484, 130.~        7.29             10.2         56.6
##  3 GB3       -0.517587, 130.~     2330.              184.         172. 
##  4 GB4       -0.517875, 130.~        8.06             11.3        192. 
##  5 GB5       -0.517319, 130.~       20.3              16.6        136. 
##  6 GB6       -0.515251, 130.~      317.               85.1        382. 
##  7 GB7       -0.515023, 130.~     1575.              150.         345. 
##  8 GB8       -0.516626, 130.~     1264.              175.         107. 
##  9 GB9       -0.514338, 130.~     1716.              152.         400. 
## 10 GB10      -0.516784, 130.~      121.               46.4        115. 
## # ... with 50 more rows, and 10 more variables: buffer_area_1000m <dbl>,
## #   tree_basal_area <dbl>, species_number <dbl>, soil_depth_mean <dbl>,
## #   leaf_litter_cover <dbl>, no_transects <dbl>, no_plots <dbl>,
## #   n_species <int>, n_genera <int>, n_family <int>

調査プロット名×種名行列（＝PCAやDCAで使える形）を作成する

countで調査プロット×植物種のクロス集計を行い、得られたデータフレームをpivot_widerで横長返還してからas.matrixで行列に変換する。その際、pivot_widerを行って生成されたNAはmutate_allとreplace_naで0に置換し、列として格納されていたplot_IDはcolumn_to_rownamesで行名に移動する。

community %>% 
  count(plot_ID, species_ID) %>% 
  pivot_wider(names_from = species_ID, values_from = n) %>% 
  mutate_all(replace_na, 0) %>% 
  column_to_rownames("plot_ID") %>% 
  as.matrix() -> species_matrix
species_matrix[1:10, 1:10]

##            species_5 species_8 species_80 species_82 species_22 species_12
## GB1_T1_ST1         1         1          1          0          0          0
## GB1_T1_ST2         2         5          0          0          0          0
## GB1_T1_ST3         0         6          0          0          0          0
## GB1_T1_ST4         0         2          1          0          0          0
## GB1_T1_ST5         0         1          0          0          0          0
## GB1_T2_ST1         0         2          0          1          0          0
## GB1_T2_ST2         0         2          0          0          1          0
## GB1_T2_ST4         0         1          0          0          0          0
## GB1_T2_ST5         0         3          0          0          0          1
## GB1_T3_ST1         0         1          0          0          0          0
##            species_2 species_63 species_26 species_43
## GB1_T1_ST1         0          0          0          0
## GB1_T1_ST2         0          0          0          0
## GB1_T1_ST3         0          0          0          0
## GB1_T1_ST4         0          0          0          0
## GB1_T1_ST5         0          0          0          0
## GB1_T2_ST1         0          0          0          0
## GB1_T2_ST2         0          0          0          0
## GB1_T2_ST4         0          0          0          0
## GB1_T2_ST5         0          0          0          0
## GB1_T3_ST1         1          4          0          0

基礎編④ggplot2を利用した作図の基礎

ggplot2を使うメリット

楽＆わかりやすい（レイヤーが明示的）＆やり直しやすい

Ex.plotしてPNGに保存

素のR：ベクトル志向

plot(iris$Sepal.Length, iris$Sepal.Width)

# PNGで保存
png("iris_scatterplot.png", width = 480, height = 480)
plot(iris$Sepal.Length, iris$Sepal.Width)
dev.off()

## png 
##   2

ggplot2：データフレーム志向

g <- ggplot() + 
          geom_point(data = iris, aes(x = Sepal.Length, y = Sepal.Width, color = Species))
plot(g)

# PNGで保存
ggsave(g, file = "iris_scatterggplot.png", width = 7, height = 7)

ggplot2の仕組みの超簡単な説明

ggplot_layer

散布図　geom_point：とりあえずプロット

ggplot2における作図はggplot関数で元レイヤを作り、geom_*関数でデータをマッピングするところから始まる。使う機会が非常に多いであろうgeom_point関数は与えられたデータを元に点をプロットする。

先に述べたように、island_propertyの60行中20行ではn_speciesの値がNAとなっている。よって、描画時にはこれらの値は描画されず、Removed 20 rows containing missing values (geom_point).という警告が出ることになる。

ggplot(data = island_property) + 
  geom_point(aes(x = island_area, y = n_species))

## Warning: Removed 20 rows containing missing values (geom_point).

g_point <- ggplot(data = island_property) + 
  geom_point(aes(x = island_area, y = n_species))

ggsave(g_point, file = "g_point.png")

## Saving 7 x 5 in image

## Warning: Removed 20 rows containing missing values (geom_point).

geom_point1

散布図　geom_point：レイヤの重ね合わせ

geom_*関数による各プロットは、レイヤとして重ね合わせることができる。

ggplot(data = island_property) + 
  geom_point(aes(x = island_area, y = n_species), color = "steelblue")

## Warning: Removed 20 rows containing missing values (geom_point).

ggplot(data = island_property) + 
  geom_point(aes(x = island_area, y = n_genera), color = "lightcoral")

## Warning: Removed 20 rows containing missing values (geom_point).

ggplot(data = island_property) + 
  geom_point(aes(x = island_area, y = n_family), color = "forestgreen")

## Warning: Removed 20 rows containing missing values (geom_point).

ggplot(data = island_property) + 
  geom_point(aes(x = island_area, y = n_species), color = "steelblue") + 
  geom_point(aes(x = island_area, y = n_genera), color = "lightcoral") + 
  geom_point(aes(x = island_area, y = n_family), color = "forestgreen")

## Warning: Removed 20 rows containing missing values (geom_point).

## Warning: Removed 20 rows containing missing values (geom_point).

## Warning: Removed 20 rows containing missing values (geom_point).

散布図　geom_point：とりあえず色分け

geom_pointに渡すaes関数内でcolor = <列名>を指定することにより、その列に含まれる値に基づいてプロット上の点の色分けを行うことができる。列に含まれる値が文字列やfactor含む離散値だった場合には色「分け」され、連続値だった場合には数値に応じたグラデーションによる描き分けがなされる。

※ここでは例示のためisland_property2を描画しているが、縦に3つ並ぶ値はすべて同じ島についてのものであることに注意。本来であればきちんとした整然データを描画することが望ましい（のだけれど、良い感じのカテゴリ値がなかったのでご容赦を）。

island_property2 <- island_property %>% 
  pivot_longer(cols = c(n_species, n_genera, n_family), names_to = "hierarchy", values_to = "value")

ggplot(data = island_property2) + 
  geom_point(aes(x = island_area, y = value, color = hierarchy)) + 
  scale_color_manual(values = c("forestgreen", "lightcoral", "steelblue"))

## Warning: Removed 60 rows containing missing values (geom_point).

geom_point3

ヒストグラム　geom_histogram

ヒストグラムの描画はgeom_histogramで行う。aesにはx軸として指定したい列を与える。階級幅は、binwidthやbreaks等のパラメータで設定する。

ggplot(data = community) + 
  geom_histogram(aes(x = tree_height_m), breaks = seq(0, 18, 1))

箱ひげ図　geom_boxplot

箱ひげ図の描画はgeom_boxplotで行う。

ggplot(data = community) + 
  geom_boxplot(aes(x = island_ID, y = tree_height_m))

棒グラフ　geom_bar

棒グラフの描画にはgeom_barを用いる。y軸の値をデータフレーム内から参照する場合にはstat = "identity"を、x軸とした値それぞれの観測数にする場合にはstat = "count"を指定する。また、積み上げ棒グラフにしたい場合にはaes内でfill = <分けたい要素が格納された列名>を指定する。割合を示したい際には、position = "fill"を指定する。

ggplot(data = community_species) + 
  geom_bar(aes(x = island_ID, fill = Family), stat = "count")

色々な微調整

scale_*関数を用いることで、x軸やy軸の目盛設定や軸ラベル設定、色分けの詳細設定などを行うことができる。

ggplot(data = island_property2) + 
  geom_point(aes(x = island_area, y = value, color = hierarchy, shape = hierarchy), alpha = 0.7) + 
  scale_x_continuous(breaks = seq(0, 12000, by = 2000)) + 
  scale_y_continuous(breaks = seq(0, 30, by = 10)) + 
  scale_color_manual(name = "Classification level", values = c("forestgreen", "lightcoral", "steelblue"), labels = c("Family", "Genera", "Species")) + 
  scale_shape_discrete(name = "Classification level", labels = c("Family", "Genera", "Species")) + 
  xlab("Island area") + 
  ylab("Number of taxa")

## Warning: Removed 60 rows containing missing values (geom_point).

other_settings

レイアウトの調整　theme

geom_関数で描画を行い、scale_関数で色分け等の微調整を行ったとしても、最終的にPNG形式等で出力した図の体裁が整うわけではない。パネル領域の背景、軸ラベルの文字、凡例の位置等の「見た目上の調整」はtheme関数内で行う必要がある。ここで注意しなければならないのが「文字や点の大きさ」であり、これらは出力画像サイズを決定してから調整することをおすすめする。

g_point_color2_theme <- ggplot(data = island_property2) + 
  geom_point(aes(x = island_area, y = value, color = hierarchy, shape = hierarchy), size = 6, alpha = 0.7) + 
  scale_x_continuous(breaks = seq(0, 12000, by = 2000)) + 
  scale_y_continuous(breaks = seq(0, 30, by = 10)) + 
  scale_color_manual(name = "Classification level", values = c("forestgreen", "lightcoral", "steelblue"), labels = c("Family", "Genera", "Species")) + 
  scale_shape_discrete(name = "Classification level", labels = c("Family", "Genera", "Species")) + 
  xlab("Island area") + 
  ylab("Number of taxa") + 
  theme(
    text = element_text(family = "sans"), 
    axis.title = element_text(size = 30, color = "black"), 
    axis.text = element_text(size = 30, color = "black"), 
    axis.line = element_line(color = "black"), 
    axis.ticks.x = element_line(color = "black"), 
    axis.ticks.length.x = unit(0.3, "cm"), 
    axis.ticks.y = element_line(color = "black"), 
    axis.ticks.length.y = unit(0.3, "cm"), 
    panel.border = element_rect(fill = NA, color = "black"), 
    plot.margin = margin(2, 2, 2, 2, "cm"), 
    legend.background = element_blank(), 
    legend.position = c(0.80, 0.15), 
    legend.title = element_text(size = 30), 
    legend.text = element_text(size = 30), 
    legend.box.background = element_blank()
  )

ggsave(g_point_color2_theme, file = "g_point_color2_theme.png", width = 14, height = 14, dpi = 200)

## Warning: Removed 60 rows containing missing values (geom_point).

g_point_color2_theme

theme1

themeは他データのプロットでも同じものを使うことができる

theme設定は出力画像サイズがほぼ変わらないのであれば、そのまま他のプロットでも使い回すことができる。これにより、複数の図のレイアウトに統一性を持たせることができる。

g_hist_theme <- ggplot(data = community) + 
  geom_histogram(aes(x = tree_height_m), breaks = seq(0, 18, 1), fill = "forestgreen") + 
  theme(
    text = element_text(family = "sans"), 
    axis.title = element_text(size = 30, color = "black"), 
    axis.text = element_text(size = 30, color = "black"), 
    axis.line = element_line(color = "black"), 
    axis.ticks.x = element_line(color = "black"), 
    axis.ticks.length.x = unit(0.3, "cm"), 
    axis.ticks.y = element_line(color = "black"), 
    axis.ticks.length.y = unit(0.3, "cm"), 
    panel.border = element_rect(fill = NA, color = "black"), 
    plot.margin = margin(2, 2, 2, 2, "cm"), 
    legend.background = element_blank(), 
    legend.position = c(0.80, 0.15), 
    legend.title = element_text(size = 30), 
    legend.text = element_text(size = 30), 
    legend.box.background = element_blank()
  )

ggsave(g_hist_theme, file = "g_hist_theme.png", width = 14, height = 14, dpi = 200)

g_hist_theme

themeはオブジェクトとして保存して呼び出すことができる&その都度上書きできる

theme設定はオブジェクトとして保存して呼び出すことができる。またすでにthemeで定義されている要素に関しても、あとからthemeで上書き再設定を行うことができる。

theme_example <- theme(
    text = element_text(family = "sans"), 
    axis.title = element_text(size = 30, color = "black"), 
    axis.text = element_text(size = 30, color = "black"), 
    axis.line = element_line(color = "black"), 
    axis.ticks.x = element_line(color = "black"), 
    axis.ticks.length.x = unit(0.3, "cm"), 
    axis.ticks.y = element_line(color = "black"), 
    axis.ticks.length.y = unit(0.3, "cm"), 
    panel.border = element_rect(fill = NA, color = "black"), 
    plot.margin = margin(2, 2, 2, 2, "cm"), 
    legend.background = element_blank(), 
    legend.position = c(0.80, 0.15), 
    legend.title = element_text(size = 30), 
    legend.text = element_text(size = 30), 
    legend.box.background = element_blank()
  )

g_hist_theme_overwrite <- ggplot(data = community) + 
  geom_histogram(aes(x = tree_height_m), breaks = seq(0, 18, 1), fill = "forestgreen") + 
  theme_example + 
  theme(
    panel.background = element_blank(), 
    panel.grid.major = element_line(color = "black", linetype = 2)
  )

ggsave(g_hist_theme_overwrite, file = "g_hist_theme_overwrite.png", width = 14, height = 14, dpi = 200)

g_hist_theme_overwrite

g_community_species <- ggplot(data = community_species) + 
  geom_bar(aes(x = island_ID, fill = Family), stat = "count") + 
  xlab("Island ID") + 
  ylab("Number of observations") + 
  theme_example + 
  theme(
    axis.text.x = element_text(angle = -45, hjust = 0), 
    legend.position = c(0.82, 0.70)
  )

ggsave(g_community_species, file = "g_community_species.png", width = 24, height = 14, dpi = 200)

g_community_species

応用：樹高とDBHを種ごとにプロットし、科で色分け

species_propertyに含まれるデータをもとに、各種の平均DBHと平均樹高をプロットし、科名で色分けする。ここではそれに加えて、geom_segmentを用いて標準誤差をもとに95%信頼区間も表示している。

g_dbh_and_height <- ggplot(data = species_perperty) + 
  geom_point(aes(x = mean_DBH, y = mean_height, color = Family), size = 4) + 
  geom_segment(aes(x = mean_DBH, xend = mean_DBH, y = mean_height - 1.96*se_height, yend = mean_height + 1.96*se_height, color = Family)) + 
  geom_segment(aes(x = mean_DBH - 1.96*se_DBH, xend = mean_DBH + 1.96*se_DBH, y = mean_height, yend = mean_height, color = Family)) + 
  scale_x_continuous(name = "Average DBH (cm)", limits = c(0, 25)) + 
  scale_y_continuous(name = "Average height (m)", limits = c(0, 15)) + 
  scale_color_discrete(name = "Family of speceis") + 
  coord_cartesian(expand = FALSE) + 
  theme_example + 
  theme(
    legend.position = c(0.80, 0.20), 
    legend.text = element_text(size = 15)
  )

ggsave(g_dbh_and_height, file = "g_dbh_and_height.png", width = 14, height = 14, dpi = 200)

## Warning: Removed 4 rows containing missing values (geom_segment).

## Warning: Removed 4 rows containing missing values (geom_segment).

g_dbh_and_height

発展編①複数グループに対して同様の処理を行う（nest→mapによる反復処理）

複雑な処理の繰り返し適用

大きなデータを解析していると複数の同形式のデータについて同様の処理を適用したいケースが出てくる
具体的には、地点ごとに生物多様性指数を算出したい、種ごとのデータに対して回帰分析を行いたい場合などが考えられる
「基礎編」で扱ったグループごとに平均値を算出するような場合はいいが、処理が複雑になった場合にも対処できる方法を知らないと、同じコードを何回もコピペしてそれぞれを訂正するという非常に危ないコーディングを行うことになる

Tidyverseを利用した繰り返し処理

Ｒの基本パッケージを用いる場合、リスト化したデータとlapplyを用いることなどが方法の一つである

iris_sp <- list(Set = iris[iris$Species == "setosa", ], ver = iris[iris$Species == "versicolor", ], vir = iris[iris$Species == "virginica", ])
correlation_sp <- lapply(iris_sp, function(x) cor.test(x$Sepal.Length, x$Sepal.Width, method = "pearson"))
correlation_sp[[1]]

## 
##  Pearson's product-moment correlation
## 
## data:  x$Sepal.Length and x$Sepal.Width
## t = 7.6807, df = 48, p-value = 6.71e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5851391 0.8460314
## sample estimates:
##       cor 
## 0.7425467

Tidyverseではtidyr::nestとpurrr::mapを用いることでより簡単に、明示的に同様の処理を行い、結果をそのままデータフレーム（Tibble）に格納することができる。

iris %>% 
　　  group_by(Species) %>% 
 　　 nest() %>% 
  　　mutate(correlation_sp = map(data, function(x) cor.test(x$Sepal.Length, x$Sepal.Width, method = "pearson")))

## # A tibble: 3 x 3
## # Groups:   Species [3]
##   Species    data              correlation_sp
##   <fct>      <list>            <list>        
## 1 setosa     <tibble [50 x 4]> <htest>       
## 2 versicolor <tibble [50 x 4]> <htest>       
## 3 virginica  <tibble [50 x 4]> <htest>

nest→mapの仕組み

nest_map

例題：島ごとに各種の観測数を算出

今回のデータで島ごとに各種の観測数を集計したい場合には次のような方法が考えられる。ここで島ごとのリストはdata2の各行の値として格納されている。

community %>% 
  group_by(island_ID) %>% 
  nest() %>% 
  mutate(
    data2 = map(data, function(data) {
                        count_of_species <- data %>% 
                                group_by(species_ID) %>% 
                                tally()
                        return(count_of_species)
                    }
)
  ) -> sp_count_per_island

sp_count_per_island

## # A tibble: 40 x 3
## # Groups:   island_ID [40]
##    island_ID data              data2            
##    <chr>     <list>            <list>           
##  1 GB1       <tibble [92 x 5]> <tibble [18 x 2]>
##  2 GB3       <tibble [84 x 5]> <tibble [17 x 2]>
##  3 GB5       <tibble [3 x 5]>  <tibble [2 x 2]> 
##  4 GB6       <tibble [45 x 5]> <tibble [8 x 2]> 
##  5 GB7       <tibble [70 x 5]> <tibble [12 x 2]>
##  6 GB8       <tibble [72 x 5]> <tibble [13 x 2]>
##  7 GB9       <tibble [91 x 5]> <tibble [18 x 2]>
##  8 GB10      <tibble [11 x 5]> <tibble [5 x 2]> 
##  9 GB11      <tibble [77 x 5]> <tibble [9 x 2]> 
## 10 GB12      <tibble [62 x 5]> <tibble [14 x 2]>
## # ... with 30 more rows

sp_count_per_island[9, 3] %>% unnest()

## Warning: `cols` is now required when using unnest().
## Please use `cols = c(data2)`

## # A tibble: 9 x 2
##   species_ID     n
##   <chr>      <int>
## 1 species_1     10
## 2 species_2      2
## 3 species_23     6
## 4 species_38     8
## 5 species_39     1
## 6 species_41     1
## 7 species_5      6
## 8 species_61     2
## 9 species_8     41

count_per_island

例題：島あたりのα多様性指数（Shannons’ H）の算出

島ごとに各種の観測数を算出してしまえば、それをもとに多様性の指標であるShannon’s Hを算出することができる。ここでは種iの観測数を島内における全種の観測数で除したものをpiとして便宜的に用いている。

※「H’」の表記もありうるが、これも便宜的に（コードを書く都合上）「H」としている。

sp_count_per_island %>% 
mutate(
    Shannon_H = map_dbl(data2, function(data2){
                            　　　percentage <- data2$n / sum(data2$n)
                            　　　H = - sum(sapply(percentage, function(x) {x * log(x)}))
                            　　　return(H)
                                }
)
  ) -> diversity_per_island

diversity_per_island

## # A tibble: 40 x 4
## # Groups:   island_ID [40]
##    island_ID data              data2             Shannon_H
##    <chr>     <list>            <list>                <dbl>
##  1 GB1       <tibble [92 x 5]> <tibble [18 x 2]>     2.40 
##  2 GB3       <tibble [84 x 5]> <tibble [17 x 2]>     2.37 
##  3 GB5       <tibble [3 x 5]>  <tibble [2 x 2]>      0.637
##  4 GB6       <tibble [45 x 5]> <tibble [8 x 2]>      1.14 
##  5 GB7       <tibble [70 x 5]> <tibble [12 x 2]>     1.91 
##  6 GB8       <tibble [72 x 5]> <tibble [13 x 2]>     2.19 
##  7 GB9       <tibble [91 x 5]> <tibble [18 x 2]>     2.25 
##  8 GB10      <tibble [11 x 5]> <tibble [5 x 2]>      1.39 
##  9 GB11      <tibble [77 x 5]> <tibble [9 x 2]>      1.54 
## 10 GB12      <tibble [62 x 5]> <tibble [14 x 2]>     2.08 
## # ... with 30 more rows

diversity_per_island

結果の描画

island_diversity <- full_join(sp_per_island, diversity_per_island, by = "island_ID") %>% full_join(island, by = "island_ID") 

g_area_diversity <- ggplot(data = island_diversity) + 
  geom_point(aes(x = island_area, y = Shannon_H), size = 6, alpha = 0.5) + 
  scale_x_continuous(name = "Island area", breaks = seq(0, 12000, by = 2000)) + 
  scale_y_continuous(name = "Shannon's H", breaks = seq(0, 3, by = 1)) + 
  theme_example
ggsave(g_area_diversity, file = "g_area_diversity.png", width = 14, height = 14, dpi = 200)

## Warning: Removed 20 rows containing missing values (geom_point).

g_distance_diversity <- ggplot(data = island_diversity) + 
  geom_point(aes(x = distance_Gam, y = Shannon_H), size = 6, alpha = 0.5) + 
  scale_x_continuous(name = "Shortest distance of each island to the nearest large landmass", breaks = seq(0, 1400, by = 200)) + 
  scale_y_continuous(name = "Shannon's H", breaks = seq(0, 3, by = 1)) + 
  theme_example
ggsave(g_distance_diversity, file = "g_distance_diversity.png", width = 14, height = 14, dpi = 200)

## Warning: Removed 20 rows containing missing values (geom_point).

g_area_diversity g_distance_diversity

発展編②ggplot2を利用した地図の描画（sfパッケージによるGISデータの描画）

ggplot2によるGISデータの描画

今回扱っているデータのように、調査地点の座標情報がわかっている場合にはデータ解析の結果を主題図のように地図上に描画することができる

Q：ただ座標をxy軸に指定したプロットじゃダメなの？
- A1：汎用性の観点からお勧めできない（ベースマップのデータはGISデータとして提供される）
- A2：座標値をそのままXY軸にしてプロットするのは妥当ではない（多くの座標情報はWGS84の緯度経度として記録されているため） ※詳しくはGISの教科書を読むか「投影座標系」とかで調べてね
Q：GISソフトウェア上で描画してはいけない？
- A：全然かまいません。ただ、GISソフトウェアは基本的にGUI操作で処理を行うので、作業の記録を残す、やり直しをしやすくするという観点からはRの方がよい（と思う）

sfパッケージでベースマップの読み込み　read_sf　データフレームのGISデータ化　st_as_f

read_sf関数を用いることで、.shpや.geojesonで保存されているGISデータをsfオブジェクトとして読み込むことができる。既存のデータフレームの座標情報を参照してそのデータをsfオブジェクトとしたい場合には、st_as_sf関数を用いる。その際には座標系をcrsに指定する必要があることに注意。

今回はベースマップとしてこちらのページで提供されているOpenStreetMapの陸域データをクリップしたものを使用した。

## install.pacages(sf)
library(sf)

## Warning: package 'sf' was built under R version 4.0.5

## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE

basemap <- read_sf("data/basemap_osm.geojson")

island_diversity_sf <- island_diversity %>% 
  separate(island_coordinates, sep = ", ", into = c("LAT", "LON")) %>% 
  mutate(LAT = str_replace(LAT, pattern = ",", replacement = "."), LON = str_replace(LON, pattern = ",", replacement = ".")) %>% 
  mutate(LAT = as.numeric(LAT), LON = as.numeric(LON)) %>% 
  st_as_sf(coords = c("LON", "LAT"), crs = "wgs84")

island_diversity_sf

## Simple feature collection with 60 features and 15 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 130.5577 ymin: -0.522189 xmax: 130.594 ymax: -0.505713
## Geodetic CRS:  WGS 84
## # A tibble: 60 x 16
##    island_ID n_species data  data2 Shannon_H island_area island_perimeter
##  * <chr>         <int> <lis> <lis>     <dbl>       <dbl>            <dbl>
##  1 GB1              18 <tib~ <tib~     2.40       4774.             303. 
##  2 GB10              5 <tib~ <tib~     1.39        121.              46.4
##  3 GB11              9 <tib~ <tib~     1.54        817.             120. 
##  4 GB12             14 <tib~ <tib~     2.08       1650.             154. 
##  5 GB13             10 <tib~ <tib~     2.12        602.              90.0
##  6 GB14             11 <tib~ <tib~     1.99        535.              90.2
##  7 GB15              9 <tib~ <tib~     1.58        381.              77.8
##  8 GB16              6 <tib~ <tib~     1.20        137.              43.5
##  9 GB17              3 <tib~ <tib~     0.509        18.4             18.1
## 10 GB18              8 <tib~ <tib~     1.78        433.              77  
## # ... with 50 more rows, and 9 more variables: distance_Gam <dbl>,
## #   buffer_area_1000m <dbl>, tree_basal_area <dbl>, species_number <dbl>,
## #   soil_depth_mean <dbl>, leaf_litter_cover <dbl>, no_transects <dbl>,
## #   no_plots <dbl>, geometry <POINT [arc_degree]>

sfオブジェクトの描画　geom_sf：とりあえずプロット

sfオブジェクトはgeom_sf関数によってggplot2レイヤーとして描画することができる。この際、座標情報はsfオブジェクトから自動的に呼び出されるため、aes内では点やポリゴンの色分けに使う列名等を指定することになる。また描画範囲はcoord_sf関数のxlimおよびylimで指定できる。

※今回の場合、island_propertyには多様性指数を算出できない（＝調査記録がない）20島のデータが含まれているため、それらの点は灰色の点となっている。

g_gis_raw <- ggplot() + 
  geom_sf(data = basemap, color = "black", fill = "darkolivegreen3") + 
  geom_sf(data = island_diversity_sf, aes(color = Shannon_H), size = 4, alpha = 0.8) + 
  scale_color_viridis_c(name = "Shannon's H", breaks = seq(0, 3, by = 1)) + 
  coord_sf(xlim = c(130.55, 130.60), ylim = c(-0.53, -0.5)) + 
  labs(caption = "Land data from OSM(https://osmdata.openstreetmap.de/data/land-polygons.html)") + 
  theme_example + 
  theme(
    legend.position = c(0.15, 0.20)
  )
ggsave(g_gis_raw, file = "g_gis_raw.png", width = 14, height = 8, dpi = 200)

g_gis_raw

投影　st_transform　と描画　geom_sf

sfオブジェクトに対してst_transform関数を適用することで座標の投影変換を行うことができる。ここでは調査地点の位置を考慮し、「UTM Zone 52S」に投影変換を行っている。投影返還を行った場合には座標単位が変わるため、coord_sfの範囲指定はそれに合わせた値を与える必要がある。ただし、描画時の軸ラベルは緯度経度表示のままとなる。

basemap_utm <- basemap %>% 
  st_transform(crs = "+proj=utm +zone=52 +south +datum=WGS84 +units=m +no_defs")

island_diversity_sf_utm <- island_diversity_sf %>% 
  st_transform(crs = "+proj=utm +zone=52 +south +datum=WGS84 +units=m +no_defs")

island_diversity_sf_utm

## Simple feature collection with 60 features and 15 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 673348.8 ymin: 9942260 xmax: 677390.6 ymax: 9944082
## CRS:           +proj=utm +zone=52 +south +datum=WGS84 +units=m +no_defs
## # A tibble: 60 x 16
##    island_ID n_species data  data2 Shannon_H island_area island_perimeter
##  * <chr>         <int> <lis> <lis>     <dbl>       <dbl>            <dbl>
##  1 GB1              18 <tib~ <tib~     2.40       4774.             303. 
##  2 GB10              5 <tib~ <tib~     1.39        121.              46.4
##  3 GB11              9 <tib~ <tib~     1.54        817.             120. 
##  4 GB12             14 <tib~ <tib~     2.08       1650.             154. 
##  5 GB13             10 <tib~ <tib~     2.12        602.              90.0
##  6 GB14             11 <tib~ <tib~     1.99        535.              90.2
##  7 GB15              9 <tib~ <tib~     1.58        381.              77.8
##  8 GB16              6 <tib~ <tib~     1.20        137.              43.5
##  9 GB17              3 <tib~ <tib~     0.509        18.4             18.1
## 10 GB18              8 <tib~ <tib~     1.78        433.              77  
## # ... with 50 more rows, and 9 more variables: distance_Gam <dbl>,
## #   buffer_area_1000m <dbl>, tree_basal_area <dbl>, species_number <dbl>,
## #   soil_depth_mean <dbl>, leaf_litter_cover <dbl>, no_transects <dbl>,
## #   no_plots <dbl>, geometry <POINT [m]>

g_gis_utm <- ggplot() + 
  geom_sf(data = basemap_utm, color = "black", fill = "darkolivegreen3") + 
  geom_sf(data = island_diversity_sf_utm, aes(color = Shannon_H), size = 4, alpha = 0.8) + 
  scale_color_viridis_c(name = "Shannon's H", breaks = seq(0, 3, by = 1)) + 
  coord_sf(xlim = c(673000, 678000), ylim = c(9942000, 9945000)) + 
  labs(caption = "Land data from OSM(https://osmdata.openstreetmap.de/data/land-polygons.html)") + 
  theme_example + 
  theme(
    legend.position = c(0.15, 0.20)
  )
ggsave(g_gis_utm, file = "g_gis_utm.png", width = 14, height = 8, dpi = 200)

g_gis_utm

最後に…

道具は道具に過ぎないから、目的に則したものを選ぼう
もしも困ったことがあったら、
- とりあえず関数のレファレンスを見よう！
- 英語で検索しよう！
- 質問サイト・ページを活用しよう！（Qiita、Slackグループ、StackOverflow）
お役立ちページ
- R for Data Science：Tidyverseの使い方についてちゃんと学ぶことができる教科書
- RStudio Cheatsheets：各パッケージの関数が図解で解説されていてめちゃくちゃ良い
連絡先：m.eco.3793（アットマーク）gmail.com

Saplings@Ecology 2022R勉強会資料（Web向け再編版）

生態学データで学ぶTidyverseの初歩

K Okano

March 19, 2022

勉強会の目的

学部生、修論生、博論生がRでよく使うデータ整形の手法について学ぶための勉強会です

テーマは「とりあえずRで何ができるかを知っておこう！」

統計処理ではなく、データ整形や集計、作図の方法についていくつか紹介します

目的（手段を増やす）に沿って研究で使いそうな処理をできるだけ使ってみます

※注意事項

勉強会の流れ

基礎編

発展編

基礎編①Rとは？Rを使う理由・利点

Rとは？

Rを使う上でとりあえず覚えておくといいこと

Rを使う理由・利点

①繰り返し＆やり直しが楽

作業面

研究面

②規格化を行いやすい

作業面

研究面

③他人に伝えるのが簡単（≒再現しやすい）

作業面

研究面

Tidyverseとは

Tidyverseを使うメリット

パイプについて

RStudioを使うメリット

RStudioとは

RStudioとは

基礎編②データ整形の基礎

準備

作業ディレクトリの設定

Tidyverseの読み込み

使うデータ：島嶼における木本植物の種多様性に関するデータ

データの読み込み read_csv

データ整形をする前に

列の選択 select

条件による行の選択 filter

ベクトルの抜き出し pull 新しい列に計算結果を格納 mutate 文字列に対する操作 str_関数

列を分割 separate 日付の定義 lubridate::ymd

観測数のカウント count

グループ化 group_by と 集計 summerise

データフレームの結合 full_join

縦長変換 pivot_longer：列名と各列の値をそれぞれ一列に

横長変換 pivot_wider：列内の要素を列名に

データフレームの書き出し write_csv

基礎編③データ整形の実践

以下のデータを得る方法を考えてみよう！

種ごとのDBHの平均と標準誤差、種ごとの樹高の平均と標準誤差と種のデータ（species、trait）を結合したデータフレームを作成する

島ごとの総出現種数、総出現属数、総出現科数と島のデータを結合したデータフレームを作成する

調査プロット名×種名行列（＝PCAやDCAで使える形）を作成する

基礎編④ggplot2を利用した作図の基礎

ggplot2を使うメリット

Ex.plotしてPNGに保存

素のR：ベクトル志向

ggplot2：データフレーム志向

ggplot2の仕組みの超簡単な説明

散布図 geom_point：とりあえずプロット

散布図 geom_point：レイヤの重ね合わせ

散布図 geom_point：とりあえず色分け

ヒストグラム geom_histogram

箱ひげ図 geom_boxplot

棒グラフ geom_bar

色々な微調整

レイアウトの調整 theme

themeは他データのプロットでも同じものを使うことができる

themeはオブジェクトとして保存して呼び出すことができる&その都度上書きできる

応用：樹高とDBHを種ごとにプロットし、科で色分け

発展編①複数グループに対して同様の処理を行う（nest→mapによる反復処理）

複雑な処理の繰り返し適用

Tidyverseを利用した繰り返し処理

nest→mapの仕組み

例題：島ごとに各種の観測数を算出

例題：島あたりのα多様性指数（Shannons’ H）の算出

結果の描画

発展編②ggplot2を利用した地図の描画（sfパッケージによるGISデータの描画）

ggplot2によるGISデータの描画

データの読み込み　read_csv

列の選択　select

条件による行の選択　filter

ベクトルの抜き出し　pull　新しい列に計算結果を格納　mutate　文字列に対する操作　str_関数

列を分割　separate　日付の定義　lubridate::ymd

観測数のカウント　count

グループ化　group_by　と　集計　summerise

データフレームの結合　full_join

縦長変換　pivot_longer：列名と各列の値をそれぞれ一列に

横長変換　pivot_wider：列内の要素を列名に

データフレームの書き出し　write_csv

散布図　geom_point：とりあえずプロット

散布図　geom_point：レイヤの重ね合わせ

散布図　geom_point：とりあえず色分け

ヒストグラム　geom_histogram

箱ひげ図　geom_boxplot

棒グラフ　geom_bar

レイアウトの調整　theme

sfパッケージでベースマップの読み込み　read_sf　データフレームのGISデータ化　st_as_f

sfオブジェクトの描画　geom_sf：とりあえずプロット

投影　st_transform　と描画　geom_sf